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Abstract 

The structure and dynamic of social network are largely 
determined by the heterogeneous interaction activity and 
social capital allocation of individuals. These features 
interplay in a non-trivial way in the formation of net¬ 
work and challenge a rigorous dynamical system theory 
of network evolution. Here we study seven real networks 
describing temporal human interactions in three differ¬ 
ent settings: scientific collaborations, Twitter mentions, 
and mobile phone calls. We find that the node’s activ¬ 
ity and social capital allocation can be described by two 
general functional forms that can be used to define a sim¬ 
ple stochastic model for social network dynamic. This 
model allows the explicit asymptotic solution of the Mas¬ 
ter Equation describing the system dynamic, and provides 
the scaling laws characterizing the time evolution of the 
social network degree distribution and individual node’s 
ego network. The analytical predictions reproduce with 
accuracy the empirical observations validating the theo¬ 
retical approach. Our results provide a rigorous dynam¬ 
ical system framework that can be extended to include 
other features of networks’ formation and to generate data 
driven predictions for the asymptotic behavior of large- 
scale social networks. 

The formation of social networks requires invest¬ 


ments in time and energy by each individual actor 
with the anticipation that collective benefits can arise 
for individuals and groups. Individuals however in¬ 
vest in developing social interactions heterogeneously 
and according to very diverse strategies. In the first 
place not all individuals are equally active in a given 
social network. Furthermore, individuals may al¬ 
locate their social capital in very diverse way, for 
instance by favoring the strengthening of a limited 
number of strong ties (bonding capital) as opposed 
to favor the exploration of weak ties opening access 
to new information and communities (bridging capi¬ 
tal) [Tj[2j[3l|4l[5l[6l[71[8]. The origins of such hetero¬ 
geneities are rooted in the trade off between compet¬ 
ing factors such as the need for close relationships [9] , 
the efforts required to keep social ties m, temporal 
and cognitive constraints HUHaiTg, and have long 
been acknowledged as key elements in the description 
of social networks’ properties [niiiiiiis], dynamical 
features [iiiiiHiiiniisEiiEaEsiEiiiiaES], and 
the the behavior of processes unfolding in social sys¬ 
tems [ii[isi[ini[i7iiini[i7iiiHi[iaEniisiiis2- How¬ 
ever, it is still lacking a general dynamical system 
framework able to relate the emerging connectivity 
pattern of social networks to the combined action of 
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social actors activity and their heterogeneity in dis¬ 
tributing resources in social capital allocation. 


Here we analyze seven time-resolved datasets de¬ 
scribing three different types of social interactions: 
scientific collaborations, Twitter mentions, and mo¬ 
bile phone calls. For all network datasets we define 
two functions statistically encoding the instantaneous 
activity of nodes and the allocation of social capital, 
respectively. The latter function is regulated by two 
parameters - system dependent- that define a simple 
reinforcement mechanism. In particular we observe 
in all datasets that the larger the number of social ties 
already activated by each node and the smaller is the 
probability of creating a new tie. We provide a thor¬ 
ough statistical characterization of the activity and 
reinforcement dynamics at play in each network and 
identify the basic parameters defining the dynamic of 
ties evolution. 


Prompted by this statistical analysis, we propose a 
dynamic network model that includes the heteroge¬ 
nous activity of nodes and the the tie formation mech¬ 
anisms. This model allows the definition of a for¬ 
mal Master Equation (ME) describing the evolution 
of the network connectivity structure that can be 
solved in the asymptotic regime (large network size 
and long time evolution). The solution of the ME 
provides the asymptotic form the degree distribution 
and the scaling relations relating degree, activity and 
and the functions characterizing the social capital al¬ 
location. The analytical solutions are capturing very 
well the empirical behavior measured in the analyzed 
datasets, connecting explicitly the evolution of social 
networks to the parameter regulating the emergence 
of heterogenous social ties. The proposed analytical 
framework is remarkably general and it can be solved 
for statistically different activity patterns. The pre¬ 
sented results have the potential to open the path to 
a general asymptotic theory of the dynamic of social 
networks by progressively integrating further social 
capital allocation strategies for the formation of so¬ 
cial ties. 


1 Results 

We analyze seven datasets containing time-stamped 
information about three different type of social inter¬ 
actions: scientific collaborations, Twitter mentions, 
and mobile phone calls. While we refer the reader 
to the Material and Methods section for the details 
of each data set, we represent all datasets as time- 
varying networks. Each node describes an individ¬ 
ual. Each time-resolved link describes a social act. 
The nature of connections is different according to 
the specific dataset. Links might represent a col¬ 
laboration resulting in a publication in a scientific 
journal, a Twitter mention, or a mobile phone call. 
We considered five scientific collaborations networks 
obtained from five different journals {PRA^ PRB, 
PRD, PRE, and PRL) of the American Physical So¬ 
ciety (APS), one Twitter mentions network (TMN), 
and one mobile phone network (MPN). 

In order to characterize the time-varying proper¬ 
ties of such networks we first measure the activity 
ai. Eormally, is defined as the fraction of inter¬ 
actions in which node i is engaged per unit of time 
with respect to all the interactions per unit time oc¬ 
curring in the network. This quantity describes the 
propensity of nodes i to be involved in social inter¬ 
actions. Empirical measurements in a wide set of 
social networks show broad distributions of activity 
[llEHlEnillSEs]. As shown in Eigurej^ [^"D], we 
confirm these observations in our datasets. In par¬ 
ticular we find that in the APS and MPN datasets 
the activity is well fitted by a truncated power law, 
while in the TMN we find a Log-Normal distribution 
(see Material and Methods and Supplementary On¬ 
line Materials for details). 

1.1 Social capital allocation 

The activity sets the clock for the activation of 
each node, however it does not provide any informa¬ 
tion on how each node invests its social capital in 
exploring new ties or reinforcing already established 
ties [34]. In order to measure the formation of new 
ties, we group nodes in classes with similar activity 
a and final degree k, so that each class b contains ac¬ 
tors with statistically equivalent characteristics (see 
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SI for details). We then measure the probability Pb{k) 
that the next social act for the nodes in the class b 
that have already contacted k nodes will result in the 
establishment of a new, k + 1-th, tie. As shown in 
Figure [E-H] Pb{k) is in general a decreasing func¬ 
tion of k. This observation resonates with previous 
research and empirical findings suggesting that our 
social interactions are bounded by cognitive and tem¬ 
poral constrains [TOl [TTJ [121 [13] . Indeed, the larger 
the number of alters in our social circle, the smaller 
the probability that the next social act will be to¬ 
wards a new tie. 

The above empirical findings suggest that the 
mechanism governing the allocation of social capital 
follow a general form that in its simplest analytical 
form can be written as: 

Pb{k)=(l + ^'^ . ( 1 ) 

In this expression, [3b modulate the tendency to ex¬ 
plore new connections, while Cb define the intrinsic 
characteristic limit of the individual to maintain mul¬ 
tiple ties. Although one could imagine more compli¬ 
cate analytical forms, we use this parsimonious ap¬ 
proach to characterize the different data sets. Inter¬ 
estingly, we find that in the five co-authorship net¬ 
works and Twitter, the exponent P is the same re¬ 
gardless of the class b. Furthermore, the values of Cb 
are typically peaked around a well defined value (see 
SI for details). More in detail, we can rescale the pro¬ 
posed functional form in each class b by defining the 
variable Xb = k/cb^ yielding 


Pb{xb)f‘ = (l+Xb) ^ (2) 

In the presence of a single exponent [3 characterizing 
the system, as shown in Figure P'K], all empirical 
curves do collapse on the reference function (l+x)“^. 
The data collapse however is not occurring in the case 
of the MPN dataset. In the latter we find a more het¬ 
erogeneous scenario in which different nodes’ classes 
are characterized by different values of [3b and c^, see 
Figure 12 L. In the Supplementary Online Material we 
provide further evidence for the evidence of a single 
or distirbuted value of p in different datasets. 


1.2 Stochastic model for the network 
dynamic 

By leveraging on the empirical evidence gathered 
here, it is possible to define a basic generative 
model of network formation based on two stochas¬ 
tic processes. Defined the network Q containing 
N nodes, at each time step a node i is active ac¬ 
cording to a probability drawn from distribution 
F{a). [23l|28l|29l|Til33]. Once active, the node i 
that has already contacted k different agents will con¬ 
tact a new, randomly chosen node with probability 
Pi{k) = (1 + k/ci)~^\ Otherwise, with probability 
1 —pi{k)^ it will interact with an already contacted 
node chosen at random. Interactions are considered 
to last one single time step. For this model it is pos¬ 
sible to write explicitly the master equation (ME) de¬ 
scribing the evolution of the probability distribution 
Pi{k^t) that a node i has degree k at time t: 


Pi{k, t -h 1) — 
P{k-l,t) 


iPi{k - 1 )+y] 

jooi k-j ^ F 


Pi{kt) 

Pi{k,t) 


ai[l-pi{k)]+Y^aj'^(l 


Pjjkj 
Piikj) 


kj 


(N-j) 


) Pj (kj , t) 
(3) 


In the above equation j ^ i and j oo i are the sum 
over the nodes already contacted and not yet con¬ 
tacted by i, respectively. Within these sums, we use 
kj as the degree of the node j. The first two terms 
on the right hand side of Eq. Q account for the cre¬ 
ation of nodes of degree k which occurs when a node 
of degree k — 1 gets active and contacts a new node, 
or when it gets in contact with a new node of previ¬ 
ous degree kj that activates and attaches to node i. 
The third and fourth terms of the r.h.s. of the equa¬ 
tion account for the conservation of nodes of degree 
k, i.e. nodes that either get active and contact one 
of their neighbors with probability a(l — p{k)) or get 
contacted by one of their neighbors. The last line of 
Eq. (§ takes into account for the case in which no 
node gets active in the current evolution time step, 
thus conserving the Pi(k,t). 
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1.3 Asymptotic theory for network with 

/3b = P 


In the case of networks characterized by a single ex¬ 
ponent (3 it is possible to consider for the ME the 
large time and large k limit, so that k can be approx¬ 
imated by a continuous variable. By neglecting the 
subleading terms of order l\t we can thus write the 
continuous asymptotic version of Eq. ^ as 


dP 

dt 


ac^ dP 

1 d^P _ 
2^ “ 


ac^ d^P a/3c^ 

I I 


P(a,/c,t)-h (4) 

dhj^P(a, h, t). 


This equation can be solved explicitly (see SI for 
details), yielding the asymptotic form: 


Pi{k, t) = ^dexp 






( 5 ) 


where ^4 is a normalization constant, C a constant 
and B{ai^Ci) a multiplicative factor of the 
term that depends on the activity and q of the 
considered agent. Its implicit expression is given in 
the SI. 

A first general result concerns the evolution in time 
of the average degree (/c(a,t)) of nodes belonging to 
a given activity class that follows the scaling laws 

(/c(a, t)) oc (at)^ . (6) 

The growth of the system is thus modulated by the 
parameter (3 that sets the strength of the reinforce¬ 
ment process in the process ruling the establishment 
of new social ties. In the limit case yd = 0 the growth 
would be linear. Indeed, the reinforcement of previ¬ 
ously activated ties would be zero and nodes would 
keep connecting randomly to other vertices, thus in¬ 
creasing indefinitely their social circle. In the oppo¬ 
site limit /3 ^ oc each node would invest is social 
capital on just one single connection, i.e. the first 
established. In the six datasets described by a sin¬ 
gle yd value, we observe the range 0.13 < yd < 0.47 
that indicates a sub-linear growth of the social sys¬ 
tem. In Eigure |S12| we find a very good agreement 
between the analytical prediction of Eq. and the 


empirical (/c(a,t)) curves, obtaining the first empir¬ 
ical validation of the modeling framework proposed 
and its ability at capturing the network formation 
dynamic. 

Eurthermore, Eq. (§ connects, at a given time t, 
the degree k and the activity a of a given node, as /c oc 
, Thus, given any specific activity distribution 
F(a), we can infer the functional form of the degree 
distribution p{k) by substituting a k^, finding: 

p{k)dk oc F{k^^^^^)k^dk. (7) 


It is important stressing that the analytical frame¬ 
work is not limited to a specific functional form of 
the activity. Indeed, with an arbitrary functional 
form of P(a), Eq. ^ gives us the possibility to pre¬ 
dict the behavior and parameters of the correspond¬ 
ing degree distribution. In Table [3]we report the 
degree distribution predicted by Eq.^^ for activities 
following a common set of heavy-tailed distributions, 
i.e. power-laws, truncated power-law, stretched ex¬ 
ponentials, and log-normal, that are usually find in 
empirical data. In Eigure |S12[ E-G] we compare the 
degree distributions p{k) predicted by Eq. (S22) with 
real data. Interestingly, also in this case the func¬ 
tional form obtained from the analytical solution of 
the model fit remarkably well the empirical evidence. 
It is important to notice that p{k) is also function 
of the parameter /3. In other words, the connectiv¬ 
ity patterns emerging from social interactions can be 
inferred knowing the propensity of individuals to be 
involved in social acts, the activity, and the strength 
of the reinforcement towards previously establish ties, 
yd. Einally it is worth remarking that Eqs. (|6, S22) 
are not affected by the distribution of q. This is an 
important result as it reduces the number of relevant 
parameters necessary to define the temporal evolu¬ 
tion of the system. 


1.4 Asymptotic theory for networks 
with distributed (3 

As we already mentioned, in the MPN dataset we 
find the evolution of social ties described by a dis¬ 
tribution of yd rather than a single value of it. This 
observation points to a more heterogeneous distribu- 
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tion of social attitudes with respect to the other six 
analyzed datasets. Arguably, such tendency might 
be driven by the different functions phone calls serve 
enabling us to communicate with relatives, friends or 
rather to companies, clients etc.. The need to intro¬ 
duce different values of (3 in the system complicates 
the model beyond analytical tractability (see SI for 
details). Nevertheless, we find that the leading term 
of the evolving average degree can be described by 
introducing a simplified model, in which the nodes of 
the system feature different values of /3 and undergo 
a simplified dynamics (see SI for further information) 
that neglects, for every node, the effects of links es¬ 
tablished by others. In these settings we can solve 
the ME and show that the minimum value of /3, /dmin, 
rules the leading term of the evolving average degree. 
In other words, we find that even in this case (/c(a, t)) 
evolves as in Eq. ® but with p substituted by /dmin- 
As shown in Eigure [ST^ D the analytical predictions 
coming from the simplified model find good agree¬ 
ment with the empirical evidence. It is interesting to 
notice that the nodes characterized by ydmin are those 
with the weak tendency to reinforce already estab¬ 
lished social ties. They are social explorers [34]. No¬ 
tably, our results, indicate that they lead the growth 
of average connectivity of the network. 

2 Discussion 

The empirical finding presented here shows clearly 
that the “cost” associated to the establishment of a 
new social tie is not constant but is function of the 
number of already activated ties, thus supporting the 
idea that social capabilities are limited by cognitive, 
temporal or other forms of constraints mi nans). 
Eraming this empirical finding in a simple stochastic 
model of network formation, we can derive a general 
asymptotic theory of the network dynamic and derive 
the general scaling laws for the behavior in time of 
the node’s degree and degree distribution. 

The model comes with some shortcomings. Indeed, 
it does not capture the modular structure or, more 
in general, correlations beyond the nearest neighbor¬ 
hood that are typical of many social networks [35] . 
In fact, individuals tend to organize their social cir¬ 


cles in tight, often hierarchical, communities. The 
model does not capture the burstiness typical of so¬ 
cial acts [ 31137 ]. We consider a simplified Poissonian 
scheme of nodes activation. A recent extension of the 
activity driven framework, without the reinforcement 
mechanism acting on social ties, has been proposed to 
account for non Poissonian node dynamics [38] . This 
is the natural starting point to generalize our model 
to bursty activities. Eurthermore, the model does 
not consider the turnover of social ties [34]. Indeed, 
in our framework once a social connection has been 
established it cannot be eliminated in favor of oth¬ 
ers. Clearly, this feature is of particular importance 
when considering social systems evolving on longer 
time scales, as the scientific journals we studied here, 
and might influence the measurement of the param¬ 
eters describing evolution of the ego-networks. 

Notwithstanding these limitations, the modeling 
framework we propose pave the way to a deeper un¬ 
derstanding of the emergence and evolution of social 
ties. The agreement between the analytical predic¬ 
tions and observed behaviors in seven real datasets, 
describing different types of social interactions, are 
encouraging steps in this direction. Einally, our re¬ 
sults are a starting point for the development of pre¬ 
dictive tools able to forecast the growth and evolution 
of social systems based not just of regression models 
or simplified toy models but on a more rigorous anal¬ 
ysis of ego-network dynamics. 

3 Materials 
3.1 Datasets 

We analyzed seven large-scale and time resolved net¬ 
works describing three different types of social inter¬ 
actions. 

• Eive networks from the APS datasets takes into 
account the co-authorship networks found in 
the Journals of the American Physical Society. 
Specifically, the PRA dataset covers the period 
from Jan. 1970 to Dec. 2006 and contains 36,880 
papers written by 34,093 authors and connected 
by 100,683 edges. The PRB dataset refers to 
the Jan. 1970 to Dec. 2007 period and con- 
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tains 104,047 papers published by 84,367 authors 
which are connected by 416,048 links. The PRD 
datasets covers the same period as the PRB one 
and it is composed by 33,376 papers, 21,202 au¬ 
thors and 60,033 edges. The PRE dataset refers 
to the Jan. 1993 to Dec. 2006 period with 24,204 
papers published by 28,188 authors connected 
by 68,029 edges. Finally, the PRL dataset con¬ 
tains ah the 66,422 papers published between 
Jan. 1960 to Dec. 2006 and written by 78,763 
authors forming 299,017 edges. 

• One network dataset describing Twitter men¬ 
tions (TMN), exchanged by users from January 
to September 2008. The network has 536,210 
nodes performing about 160M events and con¬ 
nected by 2.6M edges. 

• One Network dataset describing the mobile 
phone calls network (MPN) of 6,779,063 users of 
a single operator with about 20% market share in 
an undisclosed European country from January 
to July 2008. The datasets contains all the phone 
calls to and from company users thus including 
the calls towards or from 33,160,589 users in the 
country connected by 92,784,825 edges. 

3.2 Asymptotic solution of the ME for 
distributed values 

The solution of Eq. 0 found in Eq. 0 holds if the 
system feature a sin^e value oi p. As already dis¬ 
cussed in the MPN dataset we hnd multiple values 
of P ranging from a minimum value, /3min to a maxi¬ 
mum one /dmax- To hnd a prediction of the long time 
behavior of such a system, let us propose a simplihed 
model in which we focus on a single agent whose pa¬ 
rameters are a^, and q. In this simplihed version 
the agent can only call other nodes in the network, 
i.e. we neglect the contribution coming from the in¬ 
coming calls). In this approximation we have to solve 
a modihed version of Eq. 0 , obtained by discarding 
ah the terms containing the activity of the nodes 
j ^ i. By repeating the same procedure above, we 
get to the continuum limit that reads: 



'dPi{k,t) \d'^P{k,ty 

\k) 

dk 2 dk‘^ 


whose solution is similar to Eq. the only dif¬ 
ferences being the value of ^ and the behavior 
of the B{ai,Ci) constant (see Materials and Methods 
and the SI for details). Interestingly, even in this 

case we hnd an average degree {k{a,t)) growing ac- 

1 

cordingly to the exponent i.e. {k{a,t)) oc {at) . 

Now, let us create a reservoir of N distinct nodes of 

equal activity a and assign to each of them a diher- 

ent value of Pi drawn from an arbitrary distribution 

P{Pi). Let us also group these nodes in B classes, 

dehned so that each class i contains ah the nodes 

featuring a similar value of p Pi. If we now let 

these N nodes evolve following the simplihed model 

above, the average degree of each class i will grow as 
1 

{ki{a,t)) oc . Then, in the long time limit, the 
minimum value of pi^ i.e. /3min, will lead the growth 
of the ensemble’s average degree (see SI for further 
details), i.e. 

(A:(t)) oc t ^ + ^min . (9) 


3.3 F{a) and p{k) distributions from real 
data 

We implement the method found in [39] to determine 
the most likely functional form of both the activity 
and degree distributions. The htting procedure is as 
follows: for each functional form of the distribution 
considered (power law, log-normal, truncated power 
law and stretched exponential) we hrst determine the 
Xmin value, i.e. the lower bound to the functional 
form behavior. The Xmin value is dehned as the value 
that minimizes the Kolmgorov-Smirnov (KS) dis¬ 
tance between the analytical complementary cumu¬ 
lative distribution (CDF) and the CDF of the data. 
The latter are found for each value of Xmin by com¬ 
puting the optimal parameters of the distribution us¬ 
ing the maximum-likelihood estimator (MLE). Then, 
comparing the CDF(x > Xmin) of the data S{x) with 
the analytical one S{x)^ we compute the KS-distance 
as the maximum distance between the two CDF, i.e. 
KSdixmin) = |5'(a:;i) - P{xi)\. Once all 

the distances are computed we determine Xmin as the 
values at which the minimum distance is recorded, 
i.e. Xmin = mina, KSc^(x) (see SI and [39] for de¬ 
tails). Once we compute ah the parameters for ah 
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the functional forms analyzed we compare them with 
the likelihood ratio test IZ combined with the p-value 
that gives the statistical significance of IZ (see SI for 
details). The result of this procedure gives us the best 
candidate for the F{a) for each dataset. We find that 
a truncated power law is the best candidate for all the 
APS datasets together with the MPN one. The only 
exception is the TMN that displays a Log-Normal 
distribution of activity (see Fig. [^and SI for details). 
After we estimate the functional form and the pa¬ 
rameters of the activity distribution F(a), Eq. ( |S22| ) 
gives us the possibility to predict both the functional 
form of the degree distribution p{k) and the values 
of the parameters of such a distribution (e.g. the a 
exponent in a power-law with cutoff, see Table for 
details). The degree distribution can then be fitted 
by optimizing over the non-scale-free parameters for 
whose values we do not have an analytical or numer¬ 
ical prediction (e.g. the cut-off r in a power-law with 
cutoff). Indeed, we are missing the value of the con¬ 
stant in front of the (at)^^ term in the growth of 
the average degree (/c(a,t)) in Eq. 
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Fig. 1: (A — D) The activity distribution F(a) for PRB (^), PRL (R4, TMN (C) and the MPN (D) dataset. The 
solid lines represent the ht F(a) with the best functional form for each dataset. The latter are a truncated 
power law for the PRB, PRL and MPN case, while we find a lognormal for the TWT case (see SI and 
supplementary materials for details). In these plots we show the experimental data ranging from the lower 
bound of the fit to the 99.9% of the total amount, thus excluding the higer 0.1% of the measured activity 
values from the visible area (see materials and SI for details). (E — H) The measured pb{k) curves for selected 
nodes classes belonging to the PRB {E), PRL (T), TMN (G) and MPC {H) datasets. Each data sequence 
(different colors and markers) corresponds to a selected nodes class of the system. As one can see different 
nodes classes feature a differently behaving attachment rate function pb{k): for some nodes the probability to 
attach to a new node quickly drops to 0 at degree <10 while for some others the attachment probability is 
still > 0.1 even at very large degree {k ~ 10^). (/ — K) We rescale the attachment rate curves of all the nodes 
classes of the PRB (/), PRL (J) and TMN (K) datasets by sending k ^ Xb — kjcb and then plotting the 
Pb{xbY^^^ where {3 has the same value for every curve. For the MPC dataset (L) we show the original pb{k) 
curves belonging to a single nodes class with their ht. The resulting values of fib are shown in the legend. 
The latter are found to fall in the 1.0 < fii < 2.5 range. 
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Fig. 2: (A — D) The rescaled {k{at)) curves for selected nodes classes belonging to the PRB (A), PRL (B), TWT (C) 
and MPC (D) datasets. The time of the original data (symbols) is rescaled with the activity value t at. 
We also show the fitting curve {k{t)) oc 1 (blue solid lines) and the expected asymptotic behavior (black 
dashed lines). In the MPC case D we fit using /3 = /3min = 1-2. {E — H) The degree distribution p{k) for 
the PRB [E), PRL (T), PRA {G) and TMN {H) datasets. The predicted functional form of p(k) found in 
Eq. (S22) and Table § is shown for comparison (light blue solid lines). As in Fig. we show the data 
ranging from the lower bound of the degree distribution to the 99.9% of the data range, thus excluding from 
the plot area the higher 0 . 1 % of the degree values. 
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Fig. 3: The functional form of the activity PDF E{a) and the predicted functional form of the p{k) degree distribution 
as found in Eq. (S22), i.e. by replacing a k^~^^. This substitution fixes the scale free parameters of the 
resulting distribution, i.e. the exponent of the power-law and of the k term at the exponent in the first three 
cases, and the STD ak = 3 ^ in the Log-Normal case. The free parameters over which we fit the degree 
distribution are: (i) the cut-off r in the stretched exponential and power-law with cut-off and (ii) the 7 average 
value in the Log-Normal case. The selected PDF are, from top to bottom: power law, stretched exponential 
(Stret. Exp.), power law with cutoff (Trunc. PL) and the Log-Normal distribution. 
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Supplementary Information 


SI Data-sets 

Sl.l American Physical Society 

The APS dataset contains the five co-authorship networks of five journals of the 74nierican Physical Society, 
i.e., Physical Review T, P, P, E and Letters (P). 

The various datasets contains the data referring to all the issues of the single journals from their first issue 
up to a certain edition, specifically: 

- PRA from January 1970 to December 2006; 

- PRB and PRD from January 1970 to December 2007; 

- PRE from January 1993 to December 2006; 

- PRL from February 1960 to December 2006. 

Each dataset is composed by several files (one per month). Each file has as many lines as the number of 
papers published in that month. Einally, each line contains the IDs of the authors of the specific paper. Eor 
instance, the typical head of a file is: 

Author_000 Author_001 Author_002 #First Paper with 3 authors 

Author_003 Author_004 #Second Paper with 2 authors 


The data are cleaned so as to not take into account the papers with a single author. 

When analyzing this dataset we define the user’s activity as the number of engaged collaborations (e.g. 
an author i that publish two papers, the first with 3 co-authors and the second with a single co-author, has 
activity =4). 


SI.2 Twitter Mention Network 

The dataset of Twitter is composed by 273 daily files covering the period between January the 1®^ to 
September the 30^^2008. The dataset contains the so called fire-hose^ i.e., all the 16,329,466 citations done 
by all the 536, 210 users in the given period. The nodes in the network are connected via 2,620, 764 edges. 
Each file contains the daily events with the structure: 

Citer_ID_00 Cited_ID_00 # Event 0 

Citer_ID_01 Cited_ID_01 # Event 1 

Citer_ID_02 Cited_ID_02 # Event 2 


This dataset is not cleaned, as we have all the events that happened on the platform in the selected period. 
When analyzing this dataset we define the user’s activity as the number of citation made by i, i.e. the 
number of events actually engaged by the node i. 
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SI.3 Mobile Phone Network 

The dataset of the Mobile Phone Calls (MPC) is composed by a single file containing the 1,949,624,446 
time ordered events with 1 second resolution covering the period between January and July of 2008 for 
6, 779, 063 users of a single operator with 20% market share in an undisclosed European country. 

The dataset contains all the events from and toward users of the company (so that even the calls from 
non-company users to company users and vice-versa are taken into account). As a result, we have 33,160, 589 
nodes (of which 6, 779, 063 are users of the selected company) that are connected via 92, 784, 825 edges. 

We split the huge list of events in 98 files (each of them containing more or less the same number of events) 
for computing convenience. Each file contains events with the structure: 


Caller. 

.ID 

Called.ID 

Company_Caller 

Company_Called 

# 

Event 

0 

Caller. 

.ID 

Called.ID 

Company_Caller 

Company_Called 

# 

Event 

1 

Caller. 

.ID 

Called.ID 

Company_Caller 

Company_Called 

# 

Event 

2 


where Comp any _Caller and Company _Called are the value of the provider company of the called and caller 
nodes, respectively (e.g. the value is set to 1 if the node is a customer of our company, 0 otherwise). 

When analyzing this dataset we define the user’s activity as the number of calls done by the node, i.e. 
the number of calls actually engaged by the node i. 


S2 Data analysis 

S2.1 Activity distribution and the nodes binning 


Eor the datasets presented in Section we first evaluate, for each node i, the total number Ui of events 
engaged by the node itself. Eor instance Ui is the number of calls made by the node i in the MPC dataset 
or the number of citations done by i in the Twitter dataset. 

We then define the node activity as the ratio between the Pth node’s number of events and the total 
number of events observed in the dataset, i.e. = Uijuiot where i^tot = Uj. Thus, falls in the range 
ai G [e, 1.0) with e = mmi{ui)/utot. We then introduce and compute the activity distribution F{a). In Eig. 
SI we show the resulting activity distribution for each analyzed dataset, while in Table 0 we show the 


best candidate functional form for the F{a) distribution of each dataset. The latter is estimated using the 
methods found in [39]. 

In particular, we compare the goodness of fit on the F{a) distribution of the functional forms found in 
Table [1] of the main paper, i.e. power-law, truncated power-law, stretched exponential and log-normal 
distribution. The procedure for each dataset and each functional form reads as follows: 


- we fit the F{a) taking into account all the nodes featuring > Xmin, where Xmin is the lower bound of 
the distribution. The fit is performed using the maximum likelihood estimators (MLE) that return the 
optimal values of the parameters; 


- once the optimal parameters are found we compute the Kolmogorov-Smirnov distance {KDdixmin)) 
between the analytical and experimental complementary cumulative distribution function (CDE); 

- we then apply this procedure for different Xmin and set the amin = laina,^.^ TCSd{xmin) lower bound value 
as the one that minimizes the KSd^ 
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Fig. SI: The experimental activity distribution F{a) for (A) PRA, (B) PRB, (C) PRD, (D) PRE, (E) PRL, (E) 
TMN and (G) MPN (blue points). We also show the best candidate ht of the F{a) distribution (blue solid 
lines) featuring the functional form and parameters found in Table 0. In all the plots we show the data and 
fit ranging from the lower bound Xmin to the 99.9% of the measured data, thus excluding from the visible 
area the top 0.1% of the activity values (see Table Q for the lower bound details). 


Dataset 

Distribution 

Parameters 

KSd 

% 

£ 

TMN 

Lognormal 

<^min = 4.28e — 7, /i = —14.02, a = 1.71 

1.5e-2 

52 

—2.26e T 3 

PRA 

Truncated 

<^min = 1.90e — 5, A = 3.14e + 3, a = 1.789 

1.5e-2 

32 

-301 

PRB 

Truncated 

<^min = 6.31e — 6, A = 7.96e + 3, a = 1.638 

1.4e-2 

41 

-744 

PRD 

Truncated 

<^min = 4.54e — 5, A = 4.02e + 3, a = 1.37 

1.6e-2 

27 

-286 

PRE 

Truncated 

dmin = 4.24e — 5, A = 3.32e + 3, a = 1.92 

1.5e-2 

23 

-264 

PRL 

Truncated 

<^min = l.lOe — 5, A = 1.55e + 4, a = 1.47 

1.7e-2 

31 

-577 

MPN 

Truncated 

dmin = 2.17e — 9, A = 3.82e + 6, a = 0.448 

9.5e-3 

94 

-1.6e + 4 


Tab. 1: The candidate functional form of the activity distribution for each analyzed dataset, the evaluated parameters 
(see Table [1] in the main for the analytical expressions), the Kolmgorov-Smirnov distance KSd^ the percent 
% of nodes in the dataset that have activity ai > amin and the normalized log-likelihood C . In the parameters 
we include amin that is the value of the activity that minimizes the KS distance. This is the lower bound for 
the functional form behavior, i.e. the point at which data behave as the functional form. 
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We then repeat this procedure for all the functional forms of the F{a) and we then compare them with the 
likelihood ratio test IZ combined with the p-value that gives the statistical significance of IZ [3S1I10]. The 
result of this procedure gives us the best candidate for the F{a) for each dataset as shown in Table 0- We 
find that a truncated power law is the best candidate for all the APS datasets together with the MPN one. 
On the other hand, in the TMN we find a log-normal distribution as the best candidate for the dataset. 

Our datasets provide evidence that nodes within the same activity class (i.e. node with similar values 
of activity a^) can feature very different memory behavior. In particular agents with large activity may 
connect to very few different nodes (strong reinforcement) or establish new links at almost every step (weak 
reinforcement). For this reasons each node i of the network is naturally classified according to her activity 
Oi and her final degree ki^ i.e. the total number of different agents that have been connected to i in the 
considered time window. 

We then define a binning procedure that let us group together the similar nodes, i.e. nodes with similar 
activity and final degree. We divide the nodes in A^act activity classes so that within each activity class 
the most active node performs at most 1.5 times the events of the least active node. Then, with the same 
procedure, we further group the nodes within each activity class a according to their final degree, thus 
defining A^deg(<^) final degree classes. The nodes are therefore divided in A^deg(n) activity-degree 

classes. From now on, unless differently stated, whenever we mention the nodes’ class or bin b we will be 
referring to one of these Nij classes. 


S2.2 The reinforcement process 


To measure the reinforcement process of each system, we count all the communication events ei,{k) engaged 
by every node i of the 6-th class when it has degree ki = k. In other words, ei){k) is the total number of 
events engaged by the nodes of the 6-th class at degree k. 

Each time an event engaged by a node i of the 6-th class results in a degree increase ki = k ^ ki = k 
we increment the counter ni,{k) by 1. In other words, ni,{k) is the total number of events that the nodes 
belonging to the 6-th and featuring degree k perform toward a new node. Of course, if a node i of the 6-th 
class with degree ki = k increases its degree to ki = k 1 because it gets called by a new node, the n}y{k) 
counter is not incremented. 

The best estimate of the probability for a new node to get establish a new connection at degree k then 
reads: 


Mk) = 


nb{k) 
eb{k )’ 


(SI) 


where ni,{k) and ei,{k) are the event counters as defined above. We can give an estimate of the uncertainty 
on fb{k), by assuming that at a given degree k the events are independent (i.e. there are no correlations 
between users) and by checking that 1 <C ni,{k) <C ei,{k) so that the STD a{fi,{k)) of fb{k) reads: 


(^{fbik)) = (7b{k) = 


I fb{k){l - fb{k)) 


eb{k) 


We then fit fb{k) with the proposed reinforcement function pb{k,(3): 

-P 


(S2) 




(S3) 






S2 Data analysis 


5 



(a) 


(b) 


Fig. S2: The heat-map-like value of — ln[xb(/5)] (bottom plots). We plot the exponent j3 on the a;-axes and the 
different bins b sorted by their final degree on the j/-axes. The color-map is proportional to — ln[x^(/?)] 
representing the goodness of fit: the darker, the higher. The cyan vertical line is the value of [3opt defined 
in Eq. (S6), while the other vertical lines represent the same quantity evaluated in the three black boxes 
corresponding to different final degree intervals. (Top plots) The curve X^(/^) sis defined in Eq. (S5) (up- 
filled curve) and the same quantity for the three final degree intervals. Eor Twitter (a): a single value of 
l3opt = 0.47 fits most of the curves and only some bins b deviate from the average behavior, (b) MFC: in 
this case we observe different behaviors depending on the final degree. Thus, a single ^opt = 2.14 does not 
fit all the curves. We also show a “guide-to-the-eye” to highlight this feature (yellow dashed line). 


where c{b) is the social propensity of the 6-th bin, k is the cumulative degree and (3 is the reinforcement 
strength, that will be kept fixed for all the nodes in the system. In particular, for each class b and with a 
fixed /3, we optimize the parameter c(6), by minimizing the function xliP)- 


Kb 


XbW) = El 


[fb{k) -Pb{k,/3)Y 


k=l 


(Jb{kY 


(S4) 


where the index k runs over the Ki^ points of the 6-th bin’s curve and CFi){k) is as defined in Eq. (S2). By 
repeating this procedure for each value oi ^ ^ [0, 5.0] we find, for each class 6, a curve. 

a minimum of ^ certain 


In Fig. S2 we show the behavior of xli^)' each class 6 we find 


Popt{b) (see the horizontal lines in the heat-map-like panels of Fig. S2) 


S2 


Moreover, Fig. [S^ shows that there are two different behaviors. Specifically, in the TMN case (see Fig. 
(a)), one value of /3opt = 0.47 fits most of the curves, exception made for some outsiders: the value of /3opt(6) 
that maximizes the l/xl{P) is practically the same for all the bins. On the contrary, in the MFC case the 
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Fig. S3: Plot of the experimental curves for the (a) PRA, (b) PRD, (c) PRE, (d) PRL, (e) TMN and (f) MPC 
datasets. In the (a-e) cases the k is rescaled as k ^ k/c{b), where c{b) is the constant for the b-ih class curve 
at /3 = /3opt. The pb{k/c{b)) points are then rescaled sending In the (f) panel for MPC we 

simply plot pb{k) as a function of k with no rescaling, given that each curve features its own /3 optimal value 
/3opt as shown in the legend. 


maximum of the l/xl{P) function follows a diagonal path ranging from a larger /3opt(^) for bins with lower 
final degree to a smaller /3opt(^) for larger degree bins. In this case a single /^opt cannot fit all the curves and 
we have to consider a multi-/3 model where each class b features a different optimal value of [3^ /5opt(^)- 
In Fig. [^we present the rescaled Ph{k) curves for the PRA, PRD, PRE, PRL, TMN and MPC datasets. 
In the first five cases we show the rescaled curves obtained by substituting k k/cb and then plotting 
Pb{k) As one can see, the curves nicely collapse on the reference curve (1 In the MPC 

case we show instead the original curves, each one fitted with its own /3opt(^)- The latter parameter falls in 
the 1.2 < /3opt{b) < 3.0 interval for most of the curves as we also show in Fig. |S2[ 

To quantitatively define the /3opt parameter, let us define the total mean square deviation x^{^) 

Nb 

= '^[Xbi^)], (S5) 

6=1 

where is the total number of curves, i.e. the number of activity-degree bins b. Then, for the single 
exponent case, the function allows to define /^opt as: 

/3opt = inin;3(x^(/3)). 


(S6) 
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k [final degree] k [final degree] k [final degree] 

(a) (b) (c) 

Fig. S4: The box plot representing the distribution for different range of nodes classes b of the l3opt{b) for (a) PRB, 
(b) TMN and (c) MPN. We also show the global optimal value ^dopt (horizontal red line) as found in Eq. 
( |S6| ). The height of the box corresponds to the lower and upper quartile values of the distribution and the 
horizontal solid line corresponds to the distributions median, while the dashed lines indicates the average 
value for each range of final degree. The whiskers extend from the box to values that are within 1.5x the 
quartile range. As one can see, in both the PRB and TMN datasets the optimal values /3opt is compatible 
with the distribution found in all the nodes class ranges (we find the same result for all the other APS 
datasets analyzed). On the other hand, in the MPN the distribution of ^opt(b) lowers as the final degree 
of the class increases. The last group of nodes classes is no more compatible with the overall optimal /3opt, 
being the distribution centered around /3opt ~ IT, in agreement with our estimation of /3min = 1.2. 


In the multi-;^ case instead, we compute the different values of the exponent f3opt{b) found in the system by 
grouping the memory classes b accordingly to their final degree as shown in Fig. The optimal value of 
Popt{b) is found to be minimum for the bins featuring a large final degree, i.e. ^^min = Popt ^ 1 - 2 , which, as 
we will show in Section [S3. 2 . 3[ is the exponent driving the evolution of the network. 


To corroborate the results just outlined, we show in Fig. |^the box plot of the Popt{b) distribution for 
different groups of nodes classes b grouped by their final degree. We note that the APS and TWT datasets 
are well approximated by a single /3opt as the distribution of /3opt{b) within each sub-group of nodes is 
compatible with the global optimal value /3opt- On the other hand, in the MPN case we see that the large 
final-degree classes have their /3opt(^) distribution centered around a smaller value of /^min ^ 1-2. As already 
anticipated, this value will lead the asymptotic growth of the system as we will show in Section S3.2.3| 


As a last remark we present in Fig. (a, b) the measured distribution of the constant c{b) for the MPN 
and TMN datasets. We show the distribution for all the nodes in the network and for each activity class a, 
i.e. the group of nodes featuring similar activity. The values of this constants are distributed but peaked 
around an average value. Moreover, the distribution of the c{b) parameter within each activity closely follows 
the global one. The distribution of the social attitude c{b) then appears to be a global, activity independent 
feature of the nodes in the system. Finally, in Fig. (c) we show how the average value of the c{b) 
constants, (c) = (c( 6 )) 5 , differs from one dataset to the other varying from (c) = 0.8 in PRB to (c) = 1.7 in 
TMN and (c) = 4.6 in the MPN case, respectively. 
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(a) 


(b) 



Dataset 

Popt 

(c(6)) 

PRA 

0.20 

1.5 

PRB 

0.13 

0.8 

PRD 

0.28 

2.3 

PRE 

0.25 

2.4 

PRE 

0.15 

1.6 

TMN 

0.47 

1.7 

MPN 

1.20 

4.6 


MPN PRA PRB PRD PRE PRL TMN 


(c) 


(d) 


Fig. S5: The P{c) distribution of the constant c{h) for the (a) TMN, and (b) MPN case (solid green line). We 
also compare the global P{c) distribution with the distribution of the c{h) values found within each activity 
class (solid lines and points): we find that the distribution of the c{h) parameter is more or less activity 
independent as most of the distribution of the single activity classes follows the same functional form of the 
total distribution P(c). We then report the average value (c) of the c{h) constant for each dataset (vertical 
cyan line). The latter reads 1.71 for TMN and 4.62 for MPN. Note that in the single [3 case we evaluate c(b) 
as the values of c(b) that best fits the b-th pb(k,/S) curve fixing {3 at its optimal value (/3 = fiopt)- On the 
other hand, in the multi-^d case we evaluate c{b) as the ones that best fits the pb{k,^(b)) curve, where the 
exponent is now foxed to the 13(b) value for the memory class b, i.e. to the optimal value for the class b. (c) 
The box plot showing the global distribution of the constant c(b) in all the datasets analyzed. The height of 
the box corresponds to the lower and upper quartile values of the distribution and the horizontal solid line 
corresponds to the distributions median, while the dashed lines indicates the average value for each range of 
final degree. The whiskers extend from the box to values that are within 1.5x the quartile range, (d) In this 
table we report all the values of the reinforcement exponent ^dopt and the average reinforcement constant 
(c(6)). For the MPN case we report the /3opt = /dmin and the constant values are evaluated for each nodes 
class b using its optimal value of ^d, ydopt(^). 
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S3 The model 

S3.1 Activity driven networks with no memory 

The activity driven networks are an effective framework to describe time varying networks. The simplest 
memory-less model is defined as follows: the network consists of N nodes featuring an activity potential a^, 
i.e. the probability for a node i to get active in a certain time interval dt reads aidt. The evolution rules 
are: (i) dit each time step we start with N disconnected nodes; (ii) each node i whether gets active with 
probability aidt or does not activate with probability {1 — ai)dt. If a node gets active it calls a randomly 
selected node j in the network, thus creating an edge eij. (Hi) At the end of the time step all the created 
connection are deleted and we start again from the initial step (i). 

These evolution rules define the Master Equation (ME) for Pi{k,t), i.e. the probability that a node i of 
activity has degree k at time t, where the degree k is the number of nodes that contacted i up to time t. 
We also set, without losing generality, dt = 1. The discrete time equation for Pi{k,t) then reads: 

Pi {k^ t + 1) = 

N — k , . k , . 7-»/7 \ x—^ X—^ ^) 

ai ^ Pi{k - l,t) + ai—Pi{k,t) + Pi{k - l,t)Y^aj ^ ^ + 

jooi h 

Piik,t)'^aj'^Pj{h,t) ^^ ^ + Pi{k,t)'^aj + Pi{k,t){l - 

jooi h j^i j 

The equation is obtained in the approximation where ^ 1, so that between two consecutive times U = t 
and ti-^i = t + 1 only one site can be active. We will assume that the activity of a node i is small, i.e. 
0 < ^ 1, and we wil also consider the approximation 1 k N i.e. the integrated number of neighbors 

of a site is much larger than 1 but much smaller than the total number of agents N. The first term of the 
sum represents the probability that the site i is active and a new link is added to the system. The second 
term is the probability that the site i is active but this site connects to a site that has been already linked. 
In the third and fourth terms, the symbol denotes the sum over the sites that are not yet connected 

to i. In particular, the third term represents the probability that one of these sites is active and that it 
connects to i. The fourth term is the probability that one of these sites is active but no link between j and 
i is established. The fifth term is the probability that one of the sites already connected to i is active (being 
"^j^i the sum over the nodes already connected to i); in this case no new link is added to i. Einally, the 
last term represents the probability that at time t all the sites are not active. Eor k ^ N, the second term 
can be neglected. After some algebra we obtain the equation: 


(57) 

(58) 


Pi{k,t-\- 1) - Pi{k,t) = - {Pi{k,t) - Pi{k - 1, t)) 



given that "^^Pjih^t) = 1. Eor k we assume that = {o) i-e. the average value of the 

activity. In the limit of large time and large k we can write a continuous equation in t and k obtaining: 




— {ai + (a)) 1^- 






+ 


dk^ J 


dt 


dk 


(S9) 
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The solution of Eq. (S9) is straightforward: 


F.ik.t) = (2^(«. + + ). 


(SIO) 


In the large time limit this solution reduces to a delta function: P(a, k^t) = S{k — {a {a))t) Therefore, the 
average degree {k{a,t)) of the nodes of activity a grows as: 


(/c(a, t)) oc (a + {a))t. 


(Sll) 


as already found in mE]. Moreover the asymptotic degree distribution p{k) of a network with activity 
distribution F{a) oc a~^ is: 

p{k) (X k~^. (S12) 

S3.2 Plugging in the reinforcement process 

The model presented in Sec. |S3.1| is a basic model as it contains no correlations on an agent’s story at 
all. In particular, the probability for a node i to re-call an already contacted node is independent of the 
node degree. While simple to describe and solve analytically, this model is not realistic, as there are no 
correlations in the each agent’s history. Moreover, the probability to call an already contacted node is 
always small as k/N <C 1 (and thus the probability to call a new node remains ^ 1 even at large degree 
k). However, as shown in Sec. |S2.2 , real-world systems features a strong reinforcement process, as the 
probability Pi{k) to call a new node at degree k decreases as the degree k increases. 

For this reason we introduce an extended version of the model described in m et al. which includes a 
reinforcement function Pi{k) that measures the probability for an active node i, that has already contacted 
k different nodes in the network, to call a new node instead of an already contacted one. 


S3.2.1 The single [3 case 


As already shown in Sec. S2.2 the functional form for the reinforcement process Pi{k)^ i.e. the probability 
of adding a new link for the node i of degree k, reads: 


Pi{k) = (1 + k/ci) 


By plugging Eq. (S13) into Eq. (S8) for node i, we get; 


Pi{k,t + 1) = Pi{k - l,t) 


1) E 


jooi 


Pjjh) 

{N-h) 


Pj{h,t) 


+ 


Pi{k,t) 


i[l-pi{k)] + '^aj'^(l- 


JoOl 


Pjjh) 


N-h ^ 


Pj{h,t) 


Pi{k,t) 




(S13) 


(S14) 


where N is the number of nodes in the network, nodes not yet connected to i and 

is the sum over all the N nodes of the network. Each term of Eq. (S14) corresponds to a particular 
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event that may take place in the system, as already presented in the paper. For instance, the first term 
of the l.h.s. of Eq. (S14) takes into account the increment of the node i’s degree from k — 1 to k. This 
may happen whether because node i gets active and contacts a new node in the system with probability 
ciiPi{k — 1) or because a node j never contacted before gets active and calls exactly node i with probability 
ajPj{h)/{N — h), being h the degree of j. In the same way, the second line takes into account that node i 
does not change degree k whether because it calls an already contacted node or because the non contacted 
nodes call other nodes in the network. The last line of Eq. (S14) considers the possibility that no node in 
the network gets active. 

If we now substitute Eq. (S13) in Eq. (S14), after some algebra we get: 


Pi{k^ t + 1) — Pi{k, t) = 


dip 


{k-l + CiY 


Pi{k -l,t) - 


aic: 


{k + CiY 


Pi{k,t) 


- im.t} - p.ik - M)) g„, ^ 


(S15) 


Then, by applying the same approximations of large degree k and time t we obtain the continuous equation: 

(S16) 


dPi{k,t) _ ^cf dPi{k,t) ^ a4d^Pi{k,t) ^ 


dt 


k^ dk 2k^ dk‘^ 
ld^Pi{k,t) dPi{k,t) 

2 




dkP‘ 


dk 




where p{cj\aj) is the probability for a node j of activity aj to have reinforcement constant Cj. 


The long time asymptotic solution of Eq. ( |S16| ) is of the form: 

{k - C{ai,Ci)t^Y 


Pi{k, t) oc exp 


-yl- 


il/(l+/3) 


(S17) 


Moreover, C{a,c) is a constant depending on the activity a and the reinforcement constant c that follows 
the: 


C(a, c) 


ac^ f f ci'c'^ 


(S18) 


1 P C{a^ 

We do not have an exac t solu tion for (7(a,c), however C{a,c) for large a. 

Le t us note that Eq. (S17) can be obtained setti ng th e variable x = k — C{a)t^ and substituting it in 
Eq. (S16) and imposing that |x| ^ from Eq. (S16): 


dPi{x,t) 

dt 




C{ai,Ciy+f^t 

dPi{x, t) 


dPi{x,t) 


+ Pi(x, t) j + 


C{ai,Ci) d‘^Pi{x,t) 


dx 


dx 

f f f djdct^ 

J dajF{aj) J dcjp{cj\aj) J dy-^^-^-—^^^Pj{y,t)y. 


2(1 + 

aj/3c^ 


(S19) 


The solution of the latter equation is of the form 


Pi{x,t)^t ^<^+«exp( 


(S20) 
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thus confirming that x can be considered much smaller than 1 . 


An important consequence of equations (S17) and (S18) is that, for a system featuring a reinforcement 


strength the average degree of the nodes belonging to a class b of activity a and constant c grows as: 


(/c(a, c, t)) oc (7(a, c) • . 


(S21) 


In particular, {k{a^c^t)) oc (at)i+^ for large values of the activity a. 

As expected, the average degree grows slower than in the memoryless case {/3 = 0) where the average 
degree grows linearly in time, as found in Eq. (Sll). Moreover, the presence of a reinforcement process also 


affects the asymptotic behavior of p{k). Indeed, as already shown in the main paper, Eq. (S21) gives us 
the relation between the degree k and the activity a at a given time t, as /c oc a ^. Thus, given an activity 
distribution F(a), we can infer the functional form of the degree distribution p{k) by substituting a k^^, 
finding: 


p{k)dk oc F{k^^~^^^)k^dk. (S22) 

Specifically, by supposing a power-law activity distribution F{a) (x a~^ and considering that the degree 
distribution for a class b is described by Eq. ( |S17[ ), we obtain 

p{k) X (S23) 

where we integrated over time t and reinforcement constant c{b) and we considered the asymptotic regime 
of large time and activity. 


S3.2.2 Numerical results 

We performed numerical simulations to check the result of Section [S3. 2. 1[ We fix the following parameters: 

- N = 10^ nodes; 

- activity a G [e, 1.0] with e = 10“^, power -law distributed so that F{a) oc a~^ with u = 2.1; 

- single value of the reinforcement exponent (3 = {0.5,1.0,1.5, 2.0} and a fixed c = 1 for all the nodes; 

- T = 10^ evolution steps. 

We start with no edge in the system and we draw for each node the activity ai from the distribution F{a). 
At each step a randomly chosen node gets active with probability a^. An active node then connects with 
probability Pi{k) with a randomly chosen node which have not been yet connected to i or, with probability 
1 — pi{k), the node calls an already contacted node and no new connection is added to the system. An 
evolution step corresponds to N of these elementary steps, i.e. for each evolution step we give, on average, 
the possibility to make a call to every node in the network. 

The results are in excellent agreement with the analytical predictions. Eirst, in Eig. [^we show that the 
analysis presented in Section S2.1 correctly recovers the reinforcement exponent /dopt- Indeed the minimum 
of the xli^) vertically aligned with the value of /3 fixed in the simulations. 

Then, in Eig. S7 we present the asymptotic growth of the average degree for an activity class (i.e. a 


collection of nodes bins b featuring similar activity values) and we compare it with the analytical prediction 
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S2 


from real data. 


Fig. S6: The heat-map of — ln(x^(/3)) obtained from the simulation in the same way as in Figure ; 

As one can see the recovered ^dopt is in excellent agreement with the value used in the simulation: 1.0 in the 
(a) panel and 2.06 in the in the (b) panel. 


(/c(a, t)) (X (at) 1+/3 . In Fig. S8 
the predicted form of Eq. 


we show that the shape and the evolution of the Pi{k^t) distribution follows 


S9 


we 


The last check regards the overall degree distribution p(k) that should follow Eq. ( |S23| ). In Eig. 
compare the activity distribution F{a) oc a~^ and the degree distribution p{k). The exponent p leading the 
p{k) oc k~^ is in good agreement with the analytically predicted value /i = [(1 + P)v — /3]. 


S3.2.3 The multi-/3 case 


As shown in Section S2 a single value of the reinforcement exponent (3 is found to fit most of the Ph{^) curves 
in both the APS and TMN datasets, while for MPN a single value of P cannot fit all the Ph{^) curves for 
each activity-degree class b at once. Eor this reason we further develop the model, letting each node i to 
feature three parameters: the activity and the reinforcement constant q together with the exponent 
of the underlying reinforcement process. 

Since the model with three parameters per node is difficult to handle, we apply some approximations in 
order to get analytical insight. In particular, we work in simplified single-agent framework, where we focus 
on a single agent that can only connect to other nodes and never get called. Within this approximation the 
master equation for the node i reads: 


Pi{k, t + 1) = aip{k - l)Pi{k -l,t) + Pi{k, t) [aj(l - p{k)) + (1 - a,)]. 
The continuum limit for large degree k and time t of Eq. ds^ is: 

'dP 1 d^P] 


(S24) 


dP i'c\0 
~ ~^\k) 


dk 2 dk"^ 


(S25) 
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Fig. S7: The average degree {k{at)) for different activity classes in the {3 — 0.5 (a), 13 = 1.0 (b), /3 = 1.5 (c) and 

/3 = 2.0 (d) case. The time is rescaled with activity t at, so that all the curves collapse on a single 

1 

behavior. We also fit {k{at)) oc (t/A) (cyan solid line) and compare the simulation with the analytical 
result {k(at)) = A • t (blue dashed line). 
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Fig. S8: The probability distribution Pa{k, t) for a selected activity class a in the simulations with exponent /3 = 1.0 
(a) and /3 = 2.0 (b). We compare different evolution times (see legend) by rescaling the degree k ^ k = 
(k — (A:(a, t)))/(A;(a, on the a:-axis and the distribution itself Pa{k,t) {k{a,t)y^‘^P{a,k,t) on the 

y-axis, where {k{a, t)) is the average degree at time t for the nodes belonging to the activity class a. We also 
show the ht of the large time P{a, k, t) with a Gaussian curve (black dashed line) as predicted in Eq. ( |S17| . 


The solution for Pi{k^t) is: 


Pi{k, t) (X exp 


-A 


(k - 


21 


where the Ci now reads: 

Ci = [(1 + I3i)c^ai\ . 

Again, the average degree {ki{t)) grows as: 


{ki{t)) oc . 


(S26) 


(527) 

(528) 


The result found in Eq. ( |S28 ) holds for a single class of nodes with a given set of activity and 
reinforcement constant q and strength Pi. 

The average degree {k{aA)) for the activity class a can be computed by integrating over the different 
values of pi and q: 

{k{aA)) = j j dp'p{p', c'\a)C{a, c', Pp{t) (S29) 


where p{P,c\a) is the probability for a node of activity a to have a reinforcement exponent and constant 
equal to P and c. By assuming that the distribution of the exponent P is independent from a and c we can 
factor out the time-depend term obtaining for the activity class a: 


{k{aA)) oc J dp'p{P')t^+P 


(S30) 


where p{P) is the probability distribution of the p parameter. 
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Fig. S9: The resulting degree distribution of simulations featuring /3 = 0.5 (a), /3 = 1.0 (b), /3 = 1.5 (c) and /3 = 1.0 
(d). The analytical predictions (given F[a) oc with v — 2.1) for the scaling exponent are sown (blue 
dashed lines). 
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Let us assume that p{P) can be written as a sum of Kroenecker (^-functions, i.e.: 

Np 

pW) = 




By plugging Eq. (S31) in Eq. (S30) we find that: 


Np 




(S31) 


(S32) 




so that the minimum value of /3min leads the asymptotic behavior of the {k{a,t)) function. 

S3.2.4 Numerical results 

To investigate the multi-/d case we performed further numerical simulations considering networks with the 
following parameters: 

- N = 10^ nodes; 

- activity a G [e, 1.0] with e = 10“^, power -law distributed so that F{a) oc a~^ with z/ = 2.1; 

- (a) reinforcement exponent [3 = [0.5,1.5, 2.5] with probability [1/6,1/3,1/2] (i.e. one sixth of the nodes 
has (3 = 0.5, one third (3 = 1.5 and a half of them (3 = 2.5 regardless of their activity) and (b) 
[3 = [1.0,1.5, 2.0] with equal probability 1/3. 

- fixed c = 1 for all the nodes; 

- T = 2 • 10^ evolution steps. 

The numerical procedure is similar to the one described in Section [S3.2.2[ the difference being that we 


compute the attachment probability Pi{k) taking into account the reinforcement exponent Pi of the node 
itself. 


In Eig. |S10| we show that, in both the cases, we can recover the behavior described in Section [S2.1| for 
real data. In particular Eig. |S10[ a) (related to the P G [1,2] case) we observe a clear diagonal pattern of the 
optimal values of the exponent P{h) for the b bins that minimize the xliP)- particular P{b) varies from 
P ^ 2.0 values for lower degree nodes bins up to ^ 1.0 values for the larger final degree nodes bins. The 
figure recalls the situation of the MFC dataset presented in Eig. (b) and in the main paper. 

In Eig. Sll we show the asymptotic growth of the average degree {k{a^t)) together with the predicted 

- 1 

asymptotic behavior proportional to t i+^min . As one can see, numerical results and the suggested analytical 
solution are in very good agreement in both the cases. 

S3.3 Comparison with real data 


In Eig. |S12| we present the comparison between prediction and real data for the PRA, PRB, PRD, PRE, 
TMN and MPN datasets. The {k{a^t)) curve of each activity class is shown with the time rescaled with the 
activity of each activity class, i.e. t ^ at. In the MPC case we use as /3opt = Pmin = 1.2 (the p value found 


in the largest degree bins of Eig. S2 (b)). In all the other cases, as the /3opt fits correctly most of the curves, 
we use the /3opt value returned by our analysis. 

Einally, in Eig. |S13| we present the degree distributions, together with the predicted functional form of 
degree distribution as found in Table (1) in the main paper. 
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Fig. SIO: The heat-map like matrix of — for the simulation with /3 G [0.5,1.5, 2.5] (a) and P G [1.0,1.5, 2.0] 

(b); analogous to Figure S2 for real data. 




Fig. Sll: The average degree {k{at)) for different activity classes in the (a) 13 G [0.5,1.5.2.5] and (b) 13 G [1.0,1.5, 2.0] 

case. The time is rescaled with activity t at, so that all the curves collapse. We also plot the fit {k{at)) oc 
1 1 
(t/A) 1 +/ 3 * in the long time limit (cyan solid line) and the predicted asymptotic growth {k{at)) = AT^+^min 

(dashed line). 
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Fig. S12: The average degree {k{at)) (each data series corresponds to a different activity class) for: (a) PRA, (b) 
PRB, (c) PRD, (d) PRE, (e) TWT and (f) MPC. We compare the data for {k{a,t)) with the expected 
behavior (dashed lines) in panels (a-e) /3opt has been evaluated according Eq. (S6), while in 

the (f) case we use ^dopt = Pmin = 1-2. We also plot the power-law fit {k{a,t)) oc ^ (solid lines) 

for comparison. 
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Fig. S13: The degree distribution p{k) for: (a) PRA, (b) PRB, (c) PRD, (d) PRE, (e) PRL and (f) TMN (blue 
circles). We compare the results with the predicted behavior of Table (1) of main paper given the parameters 
of Table Q (red solid lines). We use the single value of fiopt defined by Eq. (S6) in all the cases. As in 
Eig. |S1| we show the data and fit ranging from the lower bound to the 99.9% of the measured data, thus 
excluding from the visible area the top 0.1% of the degrees values. 










